suppressWarnings(library(MendelianRandomization))
suppressWarnings(library(digest))
suppressWarnings(library(tidyr))
suppressWarnings(library(MRPRESSO))
suppressWarnings(library(rowr))
suppressWarnings(library(ggplot2))
suppressWarnings(library(dplyr))
setwd('C:/Users/marta/Desktop/IC/Thesis/Genetic Variants')
diabetes_all<-read.csv('C:/Users/marta/Desktop/IC/Thesis/Genetic Variants/Diabetes/Newrf_diabetes_instruments.csv', header = TRUE, stringsAsFactors = FALSE)
#diabetes_all<-diabetes_all[diabetes_all$F.statistic>10,] Use command to remove weak instruments
breast<-readRDS('C:/Users/marta/Desktop/IC/Thesis/Genetic Variants/Breast Cancer/breast.rds')
prostate<-readRDS('C:/Users/marta/Desktop/IC/Thesis/Genetic Variants/Prostate Cancer/prostate.rds')
breast= separate(data=breast, col=phase3_1kg_id, into=c('SNP','rest'), sep= "\\:") #Create SNP column in breast data
#Merge diabetes and breast data
index_b<-match(breast$SNP,diabetes_all$SNP..rs.id.)
diabetes_b<-cbind(diabetes_all[index_b,],breast)
#ALign effect alleles in 2 datasets
data_in_ea = diabetes_b$Effect.Allele
data_in_nea = diabetes_b$Other.Allele
ea_new = diabetes_b$a1
nea_new = diabetes_b$a0
stay=as.numeric(data_in_ea == ea_new) #stay = 1 if the same
swap=as.numeric(data_in_nea == ea_new)
#checks
table(stay)
## stay
## 0 1
## 57 54
table(swap)
## swap
## 0 1
## 54 57
stay*swap == 0
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [29] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [43] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [57] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [71] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [85] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [99] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
swap_vec = rep(1, length(swap))
swap_vec = (-2 * swap) + swap_vec
#align the beta vector in the new data
diabetes_b$bcac_onco_icogs_gwas_beta =swap_vec*as.numeric(diabetes_b$bcac_onco_icogs_gwas_beta)
diabetes_b$Effect = swap_vec * as.numeric(diabetes_b$Effect)
####################################################################################
#Merge diabetes and prostate data
index_p<-match(prostate$SNP,diabetes_all$SNP..rs.id.)
diabetes_p<-cbind(diabetes_all[index_p,],prostate)
colnames(diabetes_p)[23]<-'Effect_p'
colnames(diabetes_p)[24]<-'StdErr_p'
#ALign effect alleles in 2 datasets
data_in_ea = diabetes_p$Effect.Allele
data_in_nea = diabetes_p$Other.Allele
ea_new = diabetes_p$Allele1
nea_new = diabetes_p$Allele2
stay=as.numeric(data_in_ea == ea_new) #stay = 1 if the same
swap=as.numeric(data_in_nea == ea_new)
#checks
table(stay)
## stay
## 0
## 88
table(swap)
## swap
## 0
## 88
stay*swap == 0
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [29] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [43] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [57] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [71] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [85] TRUE TRUE TRUE TRUE
swap_vec = rep(1, length(swap))
swap_vec = (-2 * swap) + swap_vec
#align the beta vector in the new data
diabetes_p$Effect_p =swap_vec*as.numeric(diabetes_p$Effect_p)
diabetes_p$Effect = swap_vec * as.numeric(diabetes_p$Effect)
# remove NAs, we believe that they are missing at random
diabetes_b<-na.omit(diabetes_b)
diabetes_p<-na.omit(diabetes_p)
#Get rid of duplicated columns
diabetes_b <- diabetes_b[, !duplicated(colnames(diabetes_b))]
diabetes_p <- diabetes_p[, !duplicated(colnames(diabetes_p))]
#Load data of rs ids of grouped SNPs from PhenoScanner query
#Load data of rs ids of grouped SNPs of 3 molecular pathways
Lipids1<-read.csv('C:/Users/marta/Desktop/IC/Thesis/Lipids1.csv', stringsAsFactors = FALSE, header=FALSE)
Insulin1<-read.csv('C:/Users/marta/Desktop/IC/Thesis/Insulin1.csv', stringsAsFactors = FALSE, header=FALSE)
Glucose1<-read.csv('C:/Users/marta/Desktop/IC/Thesis/Glucose1.csv', stringsAsFactors = FALSE, header=FALSE)
MR breast cancer
MRInputObject_breast <- mr_input(bx = diabetes_b$Effect,
bxse = diabetes_b$StdErr,
by = diabetes_b$bcac_onco_icogs_gwas_beta,
byse = diabetes_b$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=diabetes_b$SNP)
MRAllObject_all_breast <- mr_allmethods(MRInputObject_breast, method = "all")
MRAllObject_all_breast
## Method Estimate Std Error 95% CI P-value
## Simple median -0.004 0.019 -0.041 0.033 0.835
## Weighted median 0.011 0.018 -0.024 0.045 0.539
## Penalized weighted median -0.004 0.018 -0.039 0.031 0.825
##
## IVW 0.000 0.017 -0.033 0.033 0.998
## Penalized IVW -0.007 0.013 -0.033 0.019 0.581
## Robust IVW 0.005 0.023 -0.040 0.050 0.819
## Penalized robust IVW -0.008 0.016 -0.040 0.024 0.612
##
## MR-Egger 0.029 0.035 -0.040 0.098 0.407
## (intercept) -0.003 0.003 -0.009 0.003 0.344
## Penalized MR-Egger 0.003 0.031 -0.058 0.064 0.926
## (intercept) -0.001 0.002 -0.006 0.004 0.743
## Robust MR-Egger 0.046 0.066 -0.083 0.176 0.484
## (intercept) -0.004 0.005 -0.013 0.006 0.447
## Penalized robust MR-Egger -0.005 0.041 -0.085 0.076 0.908
## (intercept) 0.000 0.003 -0.006 0.006 0.939
mr_plot(MRInputObject_breast,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast)
#MR-PRESSO breast
input = data.frame(betaY=diabetes_b$bcac_onco_icogs_gwas_beta, seY=diabetes_b$bcac_onco_icogs_gwas_se,beta_diabetes=diabetes_b$Effect, se_diabetes=diabetes_b$StdErr,row.names= NULL)
mrpresso_out = mr_presso(BetaOutcome = "betaY", BetaExposure = c("beta_diabetes"), SdOutcome = "seY", SdExposure="se_diabetes", data=input, OUTLIERtest = TRUE, DISTORTIONtest = TRUE, NbDistribution = 1000)
## Warning in mr_presso(BetaOutcome = "betaY", BetaExposure =
## c("beta_diabetes"), : Outlier test unstable. The significance threshold
## of 0.05 for the outlier test is not achievable with only 1000 to
## compute the null distribution. The current precision is <0.111. Increase
## NbDistribution.
mrpresso_out$`MR-PRESSO results`$`Global Test`$`RSSobs`
## [1] 459.3003
mrpresso_out$`MR-PRESSO results`$`Global Test`$Pvalue
## [1] "<0.001"
mrpresso_out$`Main MR results`
## Exposure MR Analysis Causal Estimate Sd T-stat
## 1 beta_diabetes Raw 4.022679e-05 0.01704120 0.002360562
## 2 beta_diabetes Outlier-corrected -7.200820e-03 0.01400442 -0.514181786
## P-value
## 1 0.9981208
## 2 0.6082263
mrpresso_out$`MR-PRESSO results`$`Distortion Test`$`Distortion Coefficient`
## beta_diabetes
## 100.5586
mrpresso_out$`MR-PRESSO results`$`Distortion Test`$Pvalue
## [1] 0.228
out_indices = mrpresso_out$`MR-PRESSO results`$`Distortion Test`$`Outliers Indices`
outliers = diabetes_b$SNP[out_indices]
# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%outliers,yes='red',no='blue'))) +
geom_point(size=1.5) +
labs(title = "", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "outliers"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 10), axis.title.x = element_text(size = 12),
axis.text.y = element_text(size = 10), axis.title.y = element_text(size = 12),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%Insulin1$V1,yes='red',no='blue'))) +
geom_point(size=2) +
labs(title = "MR analysis\n", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "Insulin SNPs"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
index<-match(diabetes_b$SNP,Insulin1$V1)
data_b_insulin<-na.omit(diabetes_b[index,])
MRInputObject_breast_insulin <- mr_input(bx = data_b_insulin$Effect,
bxse = data_b_insulin$StdErr,
by = data_b_insulin$bcac_onco_icogs_gwas_beta,
byse = data_b_insulin$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_b_insulin$SNP)
MRAllObject_all_breast_insulin <- mr_allmethods(MRInputObject_breast_insulin, method = "all")
MRAllObject_all_breast_insulin
## Method Estimate Std Error 95% CI P-value
## Simple median -0.073 0.121 -0.311 0.164 0.544
## Weighted median -0.053 0.080 -0.209 0.103 0.507
## Penalized weighted median -0.154 0.075 -0.301 -0.007 0.040
##
## IVW -0.010 0.085 -0.178 0.157 0.903
## Penalized IVW -0.108 0.074 -0.254 0.038 0.148
## Robust IVW -0.017 0.138 -0.288 0.253 0.901
## Penalized robust IVW -0.123 0.082 -0.284 0.039 0.136
##
## MR-Egger 0.012 0.141 -0.265 0.289 0.932
## (intercept) -0.002 0.009 -0.020 0.016 0.835
## Penalized MR-Egger -0.037 0.134 -0.298 0.225 0.784
## (intercept) 0.000 0.008 -0.016 0.016 0.993
## Robust MR-Egger 0.008 0.151 -0.287 0.303 0.959
## (intercept) -0.002 0.007 -0.015 0.011 0.782
## Penalized robust MR-Egger -0.039 0.121 -0.276 0.197 0.745
## (intercept) 0.000 0.006 -0.012 0.012 0.982
mr_plot(MRInputObject_breast_insulin,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_insulin)
# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%Glucose1$V1,yes='red',no='blue'))) +
geom_point(size=2) +
labs(title = "MR analysis\n", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "Glucose SNPs"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
index<-match(diabetes_b$SNP,Glucose1$V1)
data_b_glucose<-na.omit(diabetes_b[index,])
MRInputObject_breast_glucose <- mr_input(bx = data_b_glucose$Effect,
bxse = data_b_glucose$StdErr,
by = data_b_glucose$bcac_onco_icogs_gwas_beta,
byse = data_b_glucose$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_b_glucose$SNP)
MRAllObject_all_breast_glucose <- mr_allmethods(MRInputObject_breast_glucose, method = "all")
MRAllObject_all_breast_glucose
## Method Estimate Std Error 95% CI P-value
## Simple median -0.005 0.118 -0.237 0.227 0.965
## Weighted median 0.103 0.068 -0.031 0.237 0.131
## Penalized weighted median 0.209 0.068 0.075 0.343 0.002
##
## IVW 0.041 0.090 -0.135 0.217 0.650
## Penalized IVW 0.009 0.064 -0.117 0.135 0.883
## Robust IVW 0.044 0.101 -0.155 0.242 0.665
## Penalized robust IVW 0.014 0.103 -0.189 0.216 0.895
##
## MR-Egger 0.094 0.159 -0.217 0.405 0.552
## (intercept) -0.004 0.010 -0.023 0.015 0.676
## Penalized MR-Egger 0.073 0.100 -0.124 0.269 0.469
## (intercept) 0.000 0.007 -0.013 0.013 0.985
## Robust MR-Egger 0.095 0.142 -0.184 0.373 0.506
## (intercept) -0.004 0.009 -0.021 0.013 0.651
## Penalized robust MR-Egger 0.094 0.178 -0.255 0.443 0.597
## (intercept) -0.001 0.008 -0.017 0.015 0.916
mr_plot(MRInputObject_breast_glucose,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_glucose)
index<-match(diabetes_b$SNP,Lipids1$V1)
data_b_lipid<-na.omit(diabetes_b[index,])
# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%Lipids1$V1,yes='red',no='blue'))) +
geom_point(size=2) +
labs(title = "MR analysis\n", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "Lipid SNPs"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
MRInputObject_breast_lipid <- mr_input(bx = data_b_lipid$Effect,
bxse = data_b_lipid$StdErr,
by = data_b_lipid$bcac_onco_icogs_gwas_beta,
byse = data_b_lipid$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_b_lipid$SNP)
MRAllObject_all_breast_lipid <- mr_allmethods(MRInputObject_breast_lipid, method = "all")
MRAllObject_all_breast_lipid
## Method Estimate Std Error 95% CI P-value
## Simple median 0.105 0.089 -0.070 0.280 0.239
## Weighted median 0.216 0.079 0.062 0.370 0.006
## Penalized weighted median 0.226 0.082 0.065 0.388 0.006
##
## IVW 0.111 0.102 -0.089 0.311 0.276
## Penalized IVW 0.162 0.075 0.016 0.309 0.030
## Robust IVW 0.118 0.125 -0.128 0.363 0.348
## Penalized robust IVW 0.166 0.056 0.056 0.276 0.003
##
## MR-Egger -0.003 0.230 -0.453 0.447 0.989
## (intercept) 0.010 0.018 -0.026 0.047 0.571
## Penalized MR-Egger 0.093 0.222 -0.342 0.528 0.675
## (intercept) 0.006 0.017 -0.027 0.039 0.719
## Robust MR-Egger 0.000 0.243 -0.477 0.478 0.999
## (intercept) 0.010 0.020 -0.028 0.049 0.597
## Penalized robust MR-Egger 0.096 0.207 -0.309 0.502 0.641
## (intercept) 0.006 0.019 -0.030 0.042 0.745
mr_plot(MRInputObject_breast_lipid,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_lipid)
BMI=read.csv('C:/Users/marta/Desktop/IC/Thesis/Genetic variants/BMI.csv', stringsAsFactors = FALSE, header=FALSE)
index<-match(diabetes_b$SNP,BMI$V1)
data_b_bmi<-cbind(diabetes_b,BMI[index,])
data_b_bmi$V2[is.na(data_b_bmi$V2)] <- 'NO BMI'
data_nobmi<-data_b_bmi[data_b_bmi$V2=='NO BMI',]
MRInputObject_breast_nobmi <- mr_input(bx = data_nobmi$Effect,
bxse = data_nobmi$StdErr,
by = data_nobmi$bcac_onco_icogs_gwas_beta,
byse = data_nobmi$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_nobmi$SNP)
MRAllObject_all_breast_nobmi <- mr_allmethods(MRInputObject_breast_nobmi, method = "all")
MRAllObject_all_breast_nobmi
## Method Estimate Std Error 95% CI P-value
## Simple median 0.036 0.030 -0.022 0.095 0.221
## Weighted median -0.005 0.027 -0.059 0.049 0.856
## Penalized weighted median -0.004 0.028 -0.058 0.050 0.871
##
## IVW -0.012 0.022 -0.056 0.031 0.583
## Penalized IVW 0.002 0.021 -0.038 0.043 0.906
## Robust IVW -0.010 0.022 -0.052 0.033 0.649
## Penalized robust IVW 0.003 0.020 -0.036 0.041 0.898
##
## MR-Egger -0.040 0.048 -0.134 0.055 0.408
## (intercept) 0.002 0.003 -0.004 0.009 0.517
## Penalized MR-Egger -0.043 0.046 -0.133 0.048 0.358
## (intercept) 0.003 0.003 -0.004 0.009 0.374
## Robust MR-Egger -0.047 0.041 -0.127 0.034 0.254
## (intercept) 0.003 0.004 -0.004 0.011 0.400
## Penalized robust MR-Egger -0.049 0.039 -0.126 0.029 0.217
## (intercept) 0.004 0.003 -0.003 0.010 0.296
mr_plot(MRInputObject_breast_nobmi,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_nobmi)
MRInputObject_prostate <- mr_input(bx = diabetes_p$Effect,
bxse = diabetes_p$StdErr,
by = diabetes_p$Effect_p,
byse = diabetes_p$StdErr_p, outcome = 'prostate cancer', exposure='diabetes',snps=diabetes_p$SNP)
MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate, method = "all")
MRAllObject_all_prostate
## Method Estimate Std Error 95% CI P-value
## Simple median 0.011 0.024 -0.035 0.058 0.632
## Weighted median 0.001 0.022 -0.042 0.045 0.947
## Penalized weighted median -0.005 0.022 -0.049 0.039 0.819
##
## IVW -0.041 0.042 -0.124 0.042 0.333
## Penalized IVW -0.010 0.016 -0.042 0.022 0.543
## Robust IVW -0.007 0.018 -0.042 0.028 0.690
## Penalized robust IVW -0.012 0.016 -0.043 0.020 0.469
##
## MR-Egger -0.113 0.100 -0.309 0.084 0.260
## (intercept) 0.006 0.007 -0.009 0.021 0.428
## Penalized MR-Egger -0.048 0.038 -0.122 0.027 0.213
## (intercept) 0.003 0.003 -0.003 0.009 0.297
## Robust MR-Egger -0.020 0.039 -0.097 0.056 0.603
## (intercept) 0.001 0.003 -0.004 0.006 0.690
## Penalized robust MR-Egger -0.045 0.032 -0.108 0.017 0.155
## (intercept) 0.003 0.002 -0.002 0.007 0.274
mr_plot(MRInputObject_prostate,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)
#MR-PRESSO prostate
input = data.frame(betaY=diabetes_p$Effect_p,seY=diabetes_p$StdErr_p,beta_diabetes=diabetes_p$Effect,se_diabetes=diabetes_p$StdErr, row.names= NULL)
mrpresso_out = mr_presso(BetaOutcome = "betaY", BetaExposure = "beta_diabetes", SdOutcome = "seY", SdExposure="se_diabetes", data=input, OUTLIERtest = TRUE, DISTORTIONtest = TRUE, NbDistribution = 1000)
## Warning in mr_presso(BetaOutcome = "betaY", BetaExposure =
## "beta_diabetes", : Outlier test unstable. The significance threshold of
## 0.05 for the outlier test is not achievable with only 1000 to compute
## the null distribution. The current precision is <0.088. Increase
## NbDistribution.
mrpresso_out$`Main MR results`
## Exposure MR Analysis Causal Estimate Sd T-stat
## 1 beta_diabetes Raw -0.04085910 0.04223478 -0.9674276
## 2 beta_diabetes Outlier-corrected -0.01429709 0.01722820 -0.8298654
## P-value
## 1 0.3360120
## 2 0.4089954
mrpresso_out$`MR-PRESSO results`$`Distortion Test`$`Distortion Coefficient`
## beta_diabetes
## -185.7861
mrpresso_out$`MR-PRESSO results`$`Distortion Test`$Pvalue
## [1] 0.044
out_indices = mrpresso_out$`MR-PRESSO results`$`Distortion Test`$`Outliers Indices`
outliers=diabetes_p$SNP[out_indices]
ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%outliers,yes='red',no='blue'))) +
geom_point(size=1.5) +
labs(title='',x = "Genetic association with T2D", y = "Genetic association with prostate cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "outliers"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 10), axis.title.x = element_text(size = 12),
axis.text.y = element_text(size = 10), axis.title.y = element_text(size = 12),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
index<-match(diabetes_p$SNP,Insulin1$V1)
data_p_insulin<-na.omit(diabetes_p[index,])
ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%Insulin1$V1,yes='red',no='blue'))) +
geom_point(size=2) +
labs(title = "MR analysis\n", x = "T2D", y = "Prostate cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "Insulin SNPs"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
MRInputObject_prostate_insulin <- mr_input(bx = data_p_insulin$Effect,
bxse = data_p_insulin$StdErr,
by = data_p_insulin$Effect_p,
byse = data_p_insulin$StdErr_p, outcome = 'prostate cancer', exposure='T2D',snps=data_p_insulin$SNP)
MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate_insulin, method = "all")
MRAllObject_all_prostate
## Method Estimate Std Error 95% CI P-value
## Simple median 0.060 0.114 -0.164 0.285 0.598
## Weighted median 0.061 0.095 -0.125 0.247 0.520
## Penalized weighted median 0.085 0.097 -0.105 0.274 0.383
##
## IVW 0.008 0.087 -0.163 0.180 0.925
## Penalized IVW 0.062 0.082 -0.098 0.222 0.447
## Robust IVW 0.023 0.124 -0.220 0.267 0.850
## Penalized robust IVW 0.067 0.080 -0.090 0.223 0.403
##
## MR-Egger -0.001 0.178 -0.349 0.347 0.995
## (intercept) 0.001 0.011 -0.021 0.023 0.950
## Penalized MR-Egger -0.001 0.178 -0.349 0.347 0.995
## (intercept) 0.001 0.011 -0.021 0.023 0.950
## Robust MR-Egger 0.044 0.331 -0.604 0.693 0.894
## (intercept) -0.001 0.012 -0.025 0.023 0.932
## Penalized robust MR-Egger 0.044 0.331 -0.604 0.693 0.894
## (intercept) -0.001 0.012 -0.025 0.023 0.932
mr_plot(MRInputObject_prostate_insulin,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)
index<-match(diabetes_p$SNP,Glucose1$V1)
data_p_glucose<-na.omit(diabetes_p[index,])
ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%Glucose1$V1,yes='red',no='blue'))) +
geom_point(size=2) +
labs(title = "MR analysis\n", x = "T2D", y = "Prostate cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "Glucose SNPs"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
MRInputObject_prostate_glucose <- mr_input(bx = data_p_glucose$Effect,
bxse = data_p_glucose$StdErr,
by = data_p_glucose$Effect_p,
byse = data_p_glucose$StdErr_p, outcome = 'prostate cancer', exposure='diabetes',snps=data_p_glucose$SNP)
MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate_glucose, method = "all")
MRAllObject_all_prostate
## Method Estimate Std Error 95% CI P-value
## Simple median 0.163 0.110 -0.053 0.379 0.139
## Weighted median 0.094 0.088 -0.078 0.266 0.285
## Penalized weighted median 0.098 0.094 -0.086 0.281 0.297
##
## IVW 0.108 0.086 -0.061 0.276 0.210
## Penalized IVW 0.216 0.060 0.098 0.334 0.000
## Robust IVW 0.137 0.146 -0.150 0.424 0.350
## Penalized robust IVW 0.211 0.092 0.032 0.391 0.021
##
## MR-Egger 0.134 0.166 -0.191 0.459 0.418
## (intercept) -0.002 0.012 -0.025 0.021 0.849
## Penalized MR-Egger 0.335 0.109 0.121 0.549 0.002
## (intercept) -0.010 0.007 -0.023 0.004 0.164
## Robust MR-Egger 0.385 0.329 -0.260 1.030 0.242
## (intercept) -0.013 0.024 -0.060 0.035 0.599
## Penalized robust MR-Egger 0.353 0.299 -0.234 0.940 0.238
## (intercept) -0.011 0.021 -0.051 0.029 0.587
mr_plot(MRInputObject_prostate_glucose,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)
index<-match(diabetes_p$SNP,Lipids1$V1)
data_p_lipid<-na.omit(diabetes_p[index,])
ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%Lipids1$V1,yes='red',no='blue'))) +
geom_point(size=2) +
labs(title = "MR analysis\n", x = "T2D", y = "Prostate cancer", color = "Genetic variants\n") +
scale_color_manual(labels = c("Rest of SNPs", "Lipid SNPs"), values = c("blue", "red")) +
theme_bw() +
theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))
MRInputObject_prostate_lipid <- mr_input(bx = data_p_lipid$Effect,
bxse = data_p_lipid$StdErr,
by = data_p_lipid$Effect_p,
byse = data_p_lipid$StdErr_p, outcome = 'prostate cancer', exposure='diabetes',snps=data_p_lipid$SNP)
MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate_lipid, method = "all")
MRAllObject_all_prostate
## Method Estimate Std Error 95% CI P-value
## Simple median 0.090 0.100 -0.105 0.286 0.366
## Weighted median 0.089 0.089 -0.085 0.263 0.317
## Penalized weighted median 0.089 0.095 -0.098 0.276 0.351
##
## IVW 0.105 0.100 -0.090 0.301 0.291
## Penalized IVW 0.156 0.079 0.001 0.310 0.048
## Robust IVW 0.116 0.133 -0.145 0.377 0.384
## Penalized robust IVW 0.155 0.068 0.021 0.289 0.023
##
## MR-Egger 0.204 0.248 -0.282 0.689 0.410
## (intercept) -0.009 0.019 -0.046 0.029 0.659
## Penalized MR-Egger 0.542 0.143 0.262 0.823 0.000
## (intercept) -0.026 0.010 -0.046 -0.006 0.011
## Robust MR-Egger 0.571 0.037 0.498 0.645 0.000
## (intercept) -0.025 0.002 -0.029 -0.020 0.000
## Penalized robust MR-Egger 0.571 0.037 0.498 0.645 0.000
## (intercept) -0.025 0.002 -0.029 -0.020 0.000
mr_plot(MRInputObject_prostate_lipid,
error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)